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Abstract 

We report on our exploratory study for the direct evaluation of the parton dis¬ 
tribution functions from lattice QCD, based on a recently proposed new approach. 

We present encouraging results using Nf = 2 -|- 1 -|- 1 twisted mass fermions with 
a pion mass of about 370 MeV. The focus of this work is a detailed description 
of the computation, including the lattice calculation, the matching to an infinite 
momentum and the nucleon mass correction. In addition, we test the effect of 
gauge link smearing in the operator to estimate the influence of the Wilson line 
renormalization, which is yet to be done. 

1 Introduction 

Parton distribution functions (PDFs) describe the structure of hadrons by providing 
information on the momentum, angular momentum and spin of quarks and gluons in a 
hadron. Ideally, PDFs would be directly predicted by quantum chromodynamics (QCD). 
Confronted with results from deep inelastic scattering experiments, this would lead to 
a most stringent test of QCD and a deep theoretical understanding of the interaction 
between quarks and gluons. Naturally, lattice QCD methods, which can cover a broad 
energy range from the perturbative to the non-perturbative regimes, would be most 
suitable to compute the PDFs. However, such a calculation requires light-cone dynamics 
and going to short, or even zero distance on the Euclidean space-time lattice is not 
possible. 

Nevertheless, through the operator product expansion, moments of the PDFs can be 
expressed in terms of matrix elements of local operators, which are accessible to lattice 
QCD calculations. In fact, lattice QCD calculations of the PDF’s moments have been 


1 


very successfully carried out with results emerging now directly at the physical value of 
the pion mass, see the recent reviews of Refs. [D [21 a m. 

Despite the enormous activity of computing such moments in lattice QCD, it would 
still be highly desirable to have information on the PDFs themselves. A reconstruction 
of the PDFs from their moments seems unfeasible on the lattice, since higher moments 
show a very bad signal-to-noise ratio and are very hard to compute. A solution to 
this problem might be the proposal in Ref. [5], which suggests that by computing a 
parton quasidistribution function, a quantity accessible to lattice computations, contact 
to the required physical PDFs can be established through a matching procedure. Such a 
matching has already been worked out in 1-loop perturbation theory j6] and a hrst test 
of the approach has been carried out in Ref. [7] using staggered fermions. 

Here, we will use a different quark discretization on the lattice, namely twisted mass 
fermions at maximal twist [S], to conduct an exploratory study of the proposal in Ref. [S]. 
This lattice formulation of QCD has the advantage that all physical quantities scale with 
a rate of O(a^) towards the continuum limit, and it avoids the operator improvement 
necessary in other lattice QCD formulations, easing thus considerably the computa¬ 
tions. Twisted mass lattice QCD calculations for baryons have already been carried 
out successfully for the baryon spectrum Piiniiii], for form factors and moments of 
PDFs [T7[ [T31 im [T51ITB] and also for disconnected contributions to nucleon observables 

[m HH [E]. 

As stated above, our work here focuses on exploring the potential of the approach in 
Ref. j5]. To this end, we concentrate on one ensemble of maximally twisted mass fermions 
at a lattice spacing of about 0.08 fm and a pion mass of about 370 MeV. 

In our calculations, we obtain results for a boosting nucleon frame, using the three 
lowest lattice momenta, 27r/L,47r/L and Ott/L. Larger momenta show a signal-to-noise 
ratio that is too poor to extract any meaningful result. We compute the real and the 
imaginary parts of the relevant matrix elements and hnd that the imaginary part is 
very important to give an asymmetry between the quark and anti-quark distributions, 
a highly non-trivial result of our calculation. In addition, we apply different levels of 
gauge link smearing in the operator. This smearing procedure has two effects. First, 
higher smearing levels reveal the asymmetry between quark and anti-quark distribution 
much clearer. Second, different smearing levels correspond to different renormalization 
properties of the matrix elements considered. Thus, comparing results from different 
smearing levels can give a hint about the importance of renormalization, depending on 
the size of effects from smearing. We will hnally use the matching condition to relate 
the quasidistribution to the real PDF and also apply nucleon target mass correction. 

It needs to be stressed that the work presented here is only a very hrst step to 
understand the potential of the approach of Ref. j5]. It would be very important to look 
at larger momenta than used here to test that the perturbative matching works. Using a 
hypothetical mixed momentum setup (described below), we illustrate that a satisfactory 
agreement with phenomenological investigations could be obtained if larger momenta 
were available. We are planning to employ larger momenta in our next calculations by 
increasing our statistics by about an order of magnitude. 
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2 Theoretical setup 

A method to calculate quark distributions directly on a Euclidean lattice has recently 
been proposed j5]. If successful, this method can greatly improve our comprehension 
of the structure of hadrons, as well as being the hrst ab initio QCD calculation of the 
Bjorken-x dependence of the quark distributions. The key observation in this proposal 
is that from the general form of the matrix element of a twist -2 operator between a 
nucleon state with momentum P = {Pq, 0, 0, P 3 ), 


a suitable choice of the indices /xi,/i 2 , • • •,/^n makes sure that the corresponding dis¬ 
tribution is a purely spatial correlation. In Eq. ([T]), an^ are the moments of the quark 
distributions and is a symmetric rank - n tensor which can be formed with the 

target momentum P, as hrst calculated by Georgi and Politzer [20]. Let n = 2k, then 

= E(-it {g-aP-P}kAPy, (2) 


where the term {g...gP...P}kj means a symmetric sum of {2k)\/2^j\{2k — 2j)\ distinct 
products of the form p^^ 2 j+l_ ^p^^ 2 k_ Thus, setting fii = fi 2 = ■■■ = l^ 2 k = 3, 

one gets 

or 

(P|03...3|P) = 2a™(Ps)“ E 0 = 2S2.(^3)“. (4) 

with /i = /4:{P^Y M the nucleon mass. Here, we have introduced a 2 k as the matrix 

elements of the operator without subtracting the corrections in the nucleon mass. In the 
end, we want the matrix elements which can be related to the usual moments of 
the quark distributions in the Inhnite Momentum Frame (IMF). For now, we dehne 

/ +00 

x’^~^q(x,A,P3)dx, ( 5 ) 

-00 

and apply the inverse Mellin transformation to Eq. (jl]) to obtain: 

rlz 

q{x,A,P3) = J —e-^^^^P\mzh^W{z)mO)\P), (6) 

where A is the UV regulator, k^ = xP^ is the quark momentum in the ^-direction, and 
W{z) = "^ 3(2 ) jg Wilson line introduced to make the quark distribution 

gauge invariant. Eq. ([ 6 ]) is called a quasidistribution because it does not have the usual 
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properties of a quark distribution. Most notably, the momentum fraction x can be bigger 
than 1 or smaller than 0. Also, as discussed in Ref. | 6 ], the calculation of the leading UV 
divergences to the quasidistributions in perturbation theory are done keeping P 3 hxed 
while taking A —)■ cxd. This is in contrast to the case of the usual parton distributions, 
where one takes the limit P 3 —)■ 00 hrst, that is, one first goes to the IMF. The dependence 
on the UV regulator. A, will be translated, in the end, into a renormalization scale 
/i/j when relating the quasidistribution at hnite P 3 to its counterpart at inhnite P 3 . 
For now, as we still do not have a renormalization procedure for the operator and the 
coupling, we freely identify the UV regulator in the perturbative corrections in the case 
of the IMF with /x/j, the renormalization scale, while keeping it as A for the case of the 
quasidistributions. 

To relate the quasidistributions to the usual quark distributions, one uses the fact 
that the infrared r^ion of the distributions is untouched when going from a hnite to an 
inhnite momentuno. In other words, if is the usual distribution dehned though 

light-cone correlations, then one should have: 

q{x,fXR) = qbare{x)^l + ^ZF{fXR)^ + ^ j q^^\x /y, f^R^qbareiv)^ + 0{al), (7) 


q{x,A,P3) = qbare{x) \ 1 P 3 ) \ + ^ I q^^\x/y,A,P3)qbareiy)— + 


271 


a.. 


dy 


271 


' xjxc 


y 


( 8 ) 


where qbare is the bare distribution, Zp and Zp are the wave function corrections and 
and are the vertex corrections. Notice that the lower limit of integration in Eq. ([H]) 
is x/xc, where Xc ~ A/P 3 is the largest possible value of x which renders the vertex and 
wave function corrections to the quasidistributions meaningful. Opposite to the inhnite 
momentum calculation, at hnite P 3 the terms Zp and gh) ^Iq not vanish for x > 1, and 
thus this region has to be included, with the cut being made at x > 1, but below Xc- On 
the other hand, because gh)(x,/i^^) = 0 for x > 1, the integration range in Eq. ([7]) can 
be extended down to x/xc as well. 

Lattice simulations can be used to calculate the left hand side of Eq. ([S]) through 
Eq. ([ 6 ]). Ideally, one would use perturbation theory to an arbitrary order to calculate the 
right hand side of Eqs. ([7]) and ([ 8 ]) to extract the quark distribution. Currently, however, 
the self-energy and vertex corrections are known to 0{as) only and for the non-singlet 
case jB]. With this in mind, Eqs. ([7]) and dSj) can be combined to give 


g(x,A,P3) 


g(x,/iR) ^q{x,yR) ^^Zp{A,P3) - Zp^yp)^ 

+ ^ [ {q^^\x/y,A,P3)-q^^\x/y,yR))q{y,yR)— + 0{al), (9) 

jxfxc y 


and this is equivalent to Eq. (13) of Ref. j 6 ] if we consider quarks only. Notice that the 
quark number is conserved in the above expression, as long as the integrals Zp’(A, P 3 ), 

^An effective field theory approach to extract the parton distributions from the lattice observables, 
using a systematic expansion in inverse powers of the nucleon momentum, was proposed in Ref. ED- 
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listed in the Appendix, have also a cut in Xc- We dehne P 3 ) — 

Zf^Hr) and A/Pg) = gd)(^^A,P 3 ) One can include anti¬ 

quarks using the crossing relation q{x) = —q{—x), and then rewrite Eq. ([9]) as 


q{x,HR) 


q{x, A, Pg) - 
[' ZP 


X Hr A\ ^ TD\^y , rnr 
] q{y,A,P 3 )— + 0{a 

y P3 P3J y 


5 /’ 


( 10 ) 


where we have solved the system for q{x,HR). The form of Eq. fllOp that we implement 
in the actual calculations is detailed in the Appendix. 

Eq. flTOj) can be improved by calculating the corrections in M/Pg to an arbitrary 
order. As before, we write Pz)dx and use this dehnition, together 

with Eq. ([5]), to Mellin invert Eq. (jlj). After some manipulation [cf. |22]). the result is: 

9>.A) = (11) 

where f is the Nachtmann variable. The matching and the nucleon mass 

l+y 

corrections are interchangeable. 


3 Lattice calculation 

In this section, we will describe our lattice setup and our lattice computations. 

3.1 Matrix elements on the lattice 

On the lattice, the bare matrix elements h{P 3 ,z), which appear in Eq. (jS]), can be com¬ 
puted as 


h{P 3 ,z) = (P|^(z) 73 h^ 3 (A 0 )^( 0 )|P>, (12) 

with the Euclidean momentum P = (0, 0, Pg, P 4 ) and z = (0, 0, z, 0). Due to the (spatial) 
rotational symmetry on the lattice the computation can be straightforwardly applied to 
the other spatial directions. Our hnal result will then be an average over these three 
directions. 

The required matrix elements can be obtained from the ratio of suitable two- and 
three-point functions. The three-point function is constructed with the use of nucleon 
interpolating helds and a local operator; 

C'3pt(t,r,0) = (iV„(Pp)(9(r)iV«(P,0)), (13) 
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y = (y, r) 



N{P,t) 


Figure 1: Schematic picture of a possible Wick contraction of the quark fields in the three- 
point function. 


where (...) denotes the average over a sufficient number of gauge field configurations. A 
nucleon field boosted with a three-momentum can be defined via a Fourier transforma¬ 
tion of quark fields in position space: 

iV„(P, t) = Y. [d!^^{x)C-i,u%x)) , (14) 

X 

where C = *7072 and Fq,^ is a suitable parity projector. Here, we will use the parity plus 
projector F = The matrix element at vanishing momentum transfer = 0) can 

be obtained by choosing the following operator: 

0{z, t,q‘^ = 0) = Y + ^)73 1 ^ 3 ( 1 / + y)'ip{y), (15) 

y 

with y = (^, r). After Wick contracting the quark fields, the three-point function can 
be expressed in terms of quark propagators, see Fig.[T]for a schematic picture of such a 
contraction. 

We can extract the matrix element from a ratio of the above given three- and two- 
point functions: 


(F^P^(t,r,0;P) 

C'2pt(t,0;P) 


0<T<t 


-iP^ 


E 


h{P3,z), 


(16) 


where E = (Pa)^ + is the total energy of the nucleon and is the nucleon two- 

point function, which is constructed from the nucleon interpolating fields, C'^P*(f, 0; P) = 
{Na{P, t)Na{P, 0)). For the operator, we will consider the non-singlet, iso-vector quark 
combination, symbolically represented hy u — d, which avoids disconnected contributions. 

When computing the three-point function, there is a freedom on how to treat the 
propagator connecting the sink position with the operator insertion point (highlighted in 
Fig.[T]). Due to momentum projection, there is a spatial sum on both ends of the propaga¬ 
tor, which would naively require an all-to-all propagator. However, such a computation 
would need V = x T sets of inversions. 
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Here, we have tested two different methods to calculate the propagator. The first 
is the sequential method, which is exact. However, it requires the sink position and 
momentum to be fixed. As a second choice, we have used a stochastic method, where we 
use sources that contain noise on one single tinieslice (c/. Ref. |23)). The advantage 
of the stochastic method is its flexibility, allowing to freely choose the momentum at the 
sink position as well as vary the timeslice of the current insertion. 

Results from an initial test on a smaller gauge ensemble [23] indicate that both 
methods show a compatible performance and give an approximately equal error for the 
same computational effort. Thus, for the following computations the stochastic method 
will be used, since it is more flexible for studying larger momenta. 

3.2 Lattice setup 

All results shown in this work are computed on a 32^ x 64 lattice from an ETMC 
(European Twisted Mass Collaboration) production ensemble [25], with A^/ = 2 + l + l 
flavors of maximally twisted mass fermions, he. two degenerate light quarks and non¬ 
degenerate strange and charm quarks. This ensemble has a bare coupling corresponding 
to /9 = 1.95, which yields a lattice spacing of a ~ 0.082 fm m and the twisted mass 
parameter a/i = 0.0055, which corresponds to a pion mass of mps ~ 370 MeV. Our 
present statistics to compute the matrix elements is 181 gauge configurations, each with 
15 forward propagators at different source positions and two stochastic propagators, each 
propagator including both light (up and down) flavors, i. e. in total 5430 measurements. 

To examine the influence of excited states, the computation was done for two different 
source-sink separations: 8 a and 10a. From the comparison in Fig. [21 it can be seen that 
the results from both source-sink separations are visibly compatible within errors. It 
would require a significantly larger statistics to discriminate excited state effects, a task 
we want to address, however, in the future. Since here we perform an exploratory study, 
we will stick to the small separation of 8 a due to the significantly smaller noise associated 
with it. This is especially advantageous for studying larger momenta, e.g. P 3 = Gvr/L, 
which has generically a bad signal-to-noise ratio. 

3.3 Lattice results 

For our current statistics, we were able to extract matrix elements for P 3 = 27r/L,47r/L 
and Gn/L. In Ref. [7], the authors applied HYP smearing [2S] to the gauge links in the 
inserted operator. This is a lattice technique, which is used to smoothen the gauge links 
and is expected to bring the necessary renormalization factors closer to the corresponding 
tree-level value. More generally speaking, such kind of smearing will certainly influence 
the renormalization properties of the considered matrix elements. In order to obtain an 
estimate how renormalization could influence the results which will be presented here, 
we applied two and five steps of HYP smearing to the operator and compare with the 
unsmeared results in Fig. (3] 

Evidently, the effect of gauge link smearing changes the value of the matrix elements, 
for both the real and the imaginary parts. Note that the effect for the imaginary part is 
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Figure 2: We show the results for the unrenormalized matrix elements for different source-sink 
separations tg, left: P 3 = ‘I'KjL, right: P 3 = dvr/L. 


stronger than for the real part. Also, the change from zero to two steps of smearing is 
more signihcant than from two to hve steps, which indicates a saturation of the smearing 
effect. We will therefore employ a maximum of hve smearing steps in this work. We note 
in passing that a decrease of the noise like in other gluonic quantities, e.g. as in IZU, 
cannot be observed when applying smearing. A striking observation in Fig. [3] is that while 
the real part is symmetric in 2 ;, the imaginary part is highly asymmetric. This effect will 
play an important role when we discuss the quark and anti-quark distributions later on. 

Note that for a value of ^ = 0, the operator O in Eq. flT5|) can be identihed with the 
local vector current at = 0. This operator is renormalized with the vector current 
renormalization constant Zy, which, for this ensemble, is Zy = 0.625(2) [16]. After 
renormalization, the condition = 0) = 1 (c/. |2B]) should hold. Indeed, we hnd 

= 1.18(22) for P 3 = Qtt/L and Zyh““'^(0) = 0.99(3) for P 3 = Att/L while the 
value for P3 = 27r/L, Zyh““'^(0) = 0.95(1), is a bit smaller, which is probably due to 
excited state effect^. For our hnal results, we will only use data obtained for P 3 = Att/L 
and P 3 = Qtt/L. 

As can be seen in Fig. [21 when going to larger values of P 3 , the signal-to-noise ratio 
rapidly worsens. Thus, the calculation of a further, larger momentum is not possible 
with our present statistics. 

^ Using a larger source-sink separation of 10a, we find Zvh^~‘^{0) = 0.98(4), cf. Fig.j^l which is 
compatible with a value of one. We attribute the larger error to the larger source-sink separation and 
the fact that less measurements were used. 
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Figure 3: Results for the unrenormalized matrix elements with different steps of HYP 
smearing, left: real part, right: imaginary part, from top to bottom: P 3 = 
27r/L, dir/L, Gtt/L. 
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Figure 4: Comparison of results for q obtained with five, two and no steps of HYP smearing, 
Pg = 47r/L. 

4 Matching to quark distribution and nucleon mass 
corrections 

From the matrix elements P^), we calculate the quasidistributions and, after 

matching and nucleon mass corrections, the quark distributions themselves. To this end, 
we hrst apply the Fourier transformation in Eq. ([6]) to the nucleon matrix elements from 
^ = —L/2 to z = L/2, after multiplying by the vector current renormalization constant 
Zy- From this equation, it is clear that if the imaginary part of the matrix elements 
were zero, or very close to zero, there would be no difference between the positive and 
negative x regions. In other words, there would be no difference between the quark and 
antiquark distributions, as antiquarks can be interpreted as quarks in the negative x 
region, according to the crossing relation q{x) = —q{—x). 

Fig.Hlshows the complete quasidistribution for Pg = /L, after applying the Fourier 
transformation and taking the real and the imaginary parts of Pg) into account. 

An asymmetry between negative and positive x values is clearly building up, which 
is more pronounced for higher levels of gauge link smearing, emphasizing the effect 
of HYP smearing on the renormalization of these quantities Because after a proper 
renormalization the results with non-smeared and smeared gauge links have to agree 
within errors, the effect seen in Fig. H] clearly points to the fact that renormalization will 
play an important role when looking at the quark distributions obtained from lattice 
calculations in the future. 

Having the quasidistribution q{x), we can proceed to extract the physical quark dis¬ 
tribution q{x) from q{x)^ using Eq. fl2n) and then applying the nucleon mass corrections. 
To be consistent, the value of the momentum cutoff is chosen to be the same as the value 
of the lattice cutoff itself, that is, A = 1/a = 2.5 GeV. For the renormalization scale 
we make the same choice. This is a somewhat ad hoc, but plausible choice. Once a proper 
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renormalization has been carried out, the full equations for the running with fiR will be 
obtained. 

As discussed in the Appendix, the integrals also have a cut-off at Xc ~ A/P 3 , such 
that q{x > Xc, A, P 3 ) = 0. The last input we need for our calculation is the bare coupling 
constant, for which we use the value corresponding to /3 = 1.95 of our lattice calculation. 
This leads to a* = Q/{AttI3) ~ 0.245. 

We show our results in Fig. [5] for the case of 0, 2, and 5 steps of HYP smearing, 
for a nucleon with momentum P 3 = dvr/L and P 3 = Ott/L, corresponding to 0.98 GeV 
and 1.47 GeV, respectively. As anticipated, HYP smearing is essential in providing the 
required asymmetry between the quark and antiquark distributions. We note that two 
steps of smearing are already practically sufficient to show the effect of the asymmetry. 

As the nucleon momentum increases, the peak of the u{x) — d{x) distribution moves 
to smaller values of x, as it should, while for d{x) — u{x) it gets closer to 0 for most 
of the X region, but shows an increase in the small x region. This behavior is in qual¬ 
itative agreement with the behavior of the antiquark distributions as extracted from 
phenomenological analyses |291|3D1EI]. The nucleon mass corrections lead to a decrease 
of the distributions in the large x region. This is again in full accordance with our expec¬ 
tation from phenomenology and asserts that the nucleon mass corrections are essential 
to restore the energy-momentum relations, thus ensuring the partonic interpretation of 
the distributions. In addition, with increasing nucleon momentum the mass corrections 
become less and less important, as expected from Eq. ffTTj) . Finally, the slight oscillatory 
behavior in the large x region is a result of performing the Fourier transformation over 
a hnite extension only, in our case the integration is from —L/2 to +L/2. Because the 
nucleon mass corrections also make a shift of the distributions from larger to smaller 
values of x, the oscillatory behavior is more pronounced after Eq. flTT]) is applied, as the 
oscillations are more noticeable in the quasidistributions in the region a: > 1. Increasing 
the value of P3 is similar to extending the bounds of integration and thus reduces the 
oscillations. On this same line, if we had used ±1 as the limits of of integration in the 
matching, as in Eq. fIlOp and thus not taking into account the x > ^ region in the last 
two terms of Eq. (l2Tp . the oscillations would be slightly more pronounced for the case 
of P 3 = dvr/L. For P 3 = 671 /L there would be no real difference to the results presented 
in Fig.O 

Although we hud that the shape of the quark distributions resembles those of the 
phenomenological parametrizations of u{x) — d{x), with 2 or 5 steps of HYP smearing, 
we do not hud an agreement on the quantitative level. Note, however, that there is a 
clear tendency to approach the phenomenological parametrizations when P 3 is increased. 
Motivated by this observation, we made an exploratory study where we use the matrix 
elements calculated with P3 = dvr/L and P3 = bvr/L, but perform the Fourier trans¬ 
formation in Eq. ([H]), as well as the matching and the nucleon mass corrections, with 
P 3 = Stt/L. We will refer to this particular setup as the mixed momentum setup. The 
resulting distributions are shown in Fig.O 

It needs to be stressed that this exercise is, of course, only hypothetical, and using 
this setup can lead to the loss of normalization of the parton distributions, opposed to 
the results in Fig.[5l where we found a good agreement with a normalization of one. 
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Figure 5: The resulting quasidistribution q, PDF without subtracting the mass correction q, 
and hnal PDF q^^\ left: P 3 = dvr/L, right: P 3 = Gvr/L, from top to bottom: 0 
steps, 2 steps, 5 steps of HYP smearing, negative region: q{x) = —q{—x), comparison 
with phenomenological u{x) — d{x) curves at = 6.25MeV^ (MSTW [29], CJ12 
[30], ABMll |3l]). 
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Figure 6: Results from an hypothetical mixed momentum analysis using different values of 
the momentum in the computation of the lattice matrix element (left: P 3 = dvr/L, 
right: P 3 = Gvr/L) than in the Fourier transformation, matching and mass correc¬ 
tions (P 3 = Stt/L) with 5 steps of HYP smearing. 


Nevertheless, the agreement with the phenomenological parametrizations of the distri¬ 
butions at the intermediate and large x regions is really encouraging. This indicates 
that by employing an only moderately larger value of P 3 than the ones used here, we 
could obtain even a quantitative agreement to the parametrizations in certain regions 
of X . This concerns in particular the large x region, where increasing values of P3 tend 
to bring the resulting distribution down. In the small (and positive) x region, on the 
other hand, it seems that increasing the nucleon momentum is not sufficient to produce 
a rise of the distribution. This may be related to the fact that there is a limitation in 
the present calculation in the small x region due to the presence of the infrared, 1/L, 
and ultra-violet, 1/a cut-off regulators on a finite lattice. Thus, this limitation will be 
overcome when larger lattices and smaller values of the lattice spacing become avail¬ 
able. Furthermore, we stress that the here obtained results are at only one, non-physical 
value of the pion mass and the shape of distribution might as well depend on the quark 
mass. In any case, a more definite statement can only be made after we have access 
to the matrix elements for P 3 = Stt/L. This is not possible with our present statistics. 
However, we are in the process of generating a substantially higher statistics. This will 
allow us to extrapolate the data for P3 = 27 r/L, P3 = 47r/L and P3 = bvr/L to obtain 
the quasidistribution at P 3 = Stt/L. Although we do not expect a big difference to the 
situation of the hypothetical mixed setup shown in Fig. [ 6 ], a full analysis with real data 
is, of course, mandatory and will be presented in a forthcoming work. 


5 Conclusions and outlook 

In this work, we have presented our first exploratory study of the approach developed 
in Ref. [5] for the calculation of the x dependence of quark distributions directly on the 


13 



















lattice, employing the twisted mass formulation of lattice QCD. The study presented 
here, together with the work of Ref. [7|, constitute the first two attempts to implement 
the approach of Ref. [5] in realistic lattice QCD calculations. Our results, represented in 
Fig.O are comparable to those of Fig. 2 of [7] and we see that the two calculations give 
similar results for the case of 2 steps of HYP smearing. Yet it seems that in our case, for 
Fs = hvr/L, the shift of the peak of the quark distributions towards the small x region 
is bigger. Also, as we increase the number of steps of HYP smearing, the position of the 
peaks is basically unchanged, but they are more pronounced. On the other hand, the 
large x region tends to be smaller. It is very reassuring to see that both effects bring 
our results closer to the phenomenological parametrizations. Beyond these effects, HYP 
smearing is fundamental to generate a sizeable (asymmetric in z) imaginary part in the 
matrix elements. This result generates automatically an asymmetry between the quark 
and the antiquark distributions, a highly non-trivial result. 

The outcome of our ab initio lattice QCD calculation with a small and positive 
d{x) — u{x) is in a very good qualitative agreement with phenomenological parametriza¬ 
tions. In a hypothetical exercise where we use a larger momentum of F 3 = Svr/L in the 
Fourier transformation than we actually have in our lattice QCD calculation, we observe 
a better qualitative behavior as compared to what is expected phenomenologically, as 
is shown in Fig.El Moreover, it is clear from both Figs.|5] and E] that increasing the mo¬ 
mentum implies only marginal corrections to the quasidistributions, the corrections for 
the case P 3 = Gtt/L being restricted from intermediate to small x regions only. 

In summary, we have presented our first effort to explore the potential to calculate 
quark distributions directly within the lattice QCD formulation. Although there are 
clearly shortcomings, such as not being able to reach large momenta and the lack of 
renormalization, our results are promising. In particular, our study of the quark distri¬ 
bution in the mixed momentum setup indicates that only moderately larger momenta 
than used here may be sufficient to reach a quantitative agreement with phenomeno¬ 
logical parametrizations in the large x region. We are presently increasing our statistics 
signihcantly, which will allow us to obtain data with such larger momenta. In addition, 
we are testing different approaches to perform the necessary renormalization of the ma¬ 
trix elements entering the calculation of the quasidistributions. Finally, applying this 
method to gluon conhgurations generated directly at the physical value of the pion mass 
may open the exciting possibility to address quark distributions and therefore unravel 
the structure of hadron from first principle QCD calculations. 
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Appendix 


The wave function and vertex corrections in Eq. (IlOp were calculated in Ref. |S]. The 
vertex corrections are given by: 
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for < 0 . The wave function corrections are given by: 
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for .^ < 0. 

In the actual calculation, we make a change of variables in the integral term contain¬ 
ing q{y, A, P3) of Eq. flT0|) . and also set the threshold above which the quasidistribution is 
zero. We call this value Xc, which is of order of A/P3. When we inverted Eq. ([ 9 ]), we kept 
the limits of integration from -1 to +1, which is the region where the quark distributions 
are defined and where factorization holds. In practice, we will integrate from —Xc to +Xc, 
the reason being that q{x > 1) 7^ 0 and contributions from this region should be taken 
into account. As we increase the value of P3, however, the closer we get to the physical 
distribution and as a result q{x > 1 ) ~ 0 . We also break the integral containing q into 
two terms, with the limits from —Xc to — \x\/xc and from + \x\/xc to +Xc- We then make 
a change of variables, ^ = x/y, and Eq. flTOl) is rewritten as: 
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q(,x,fiR) 



The integrals contain both single and double poles at .^ = 1. It can be shown that 
the single pole terms cancel between Eqs. ffT7|) - ffT8|) and (l^ - fl22|) . e.g. the single pole 
in the third term on the l.h.s of Eq. 0181) is cancelled by the third term on the hh.s of 
Eq. fl22D . The double poles are first reduced to a single pole by a similar cancellation 
when combining the vertex and wave function corrections, as in the single pole case, and 
the remaining pole is taken care of by using the Cauchy’s principal value prescription. 
The remaining expression is finite, with the exception that the integral of is 

divergent as ^ ^ ±oo. The divergent term is: 


q{x)^\n{xl-l), (25) 

where we have set Xc as the upper and lower limit of the integrals of fl2Tl) and fl2^ . 
respectively. The same limits of integration, both when integrating and SZ^^'^ are 
necessary in order to maintain the quark number conservation. Notice that this log 
divergent term is the usual UV divergence present in the wave function. 
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